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1 Introduction 


We have developed a new method to study nonstationary signals and systems. Modern engineering 
has developed the concept and methods generally called “system theory” or “input-output rela¬ 
tions” and the associated concepts. This has become textbook material. We believe we generalized 
the concept of “input-output relations” by developing a new method to study systems and channels 
that applies to both the time-variant and time-invariant case. We have shown that a remarkable 
simplification occurs, both conceptually and technically, when systems are formulated in phase 
space. The phase space can be time-frequency for the usual formulation of channels/systems, but 
can also be position-wave number, for systems such as pulse propagation and array processing. 
Our method has led to practical solutions to many problems, both deterministic and random. The 
fundamental idea of our approach can be best conveyed by a contrast with the traditional method. 
The standard textbook approach which has been the foundation of input-output relations is the 
concept of a “system function”, or Green’s function, or impulse response function. If we have an 
input time function, f(t) : which passes through a system characterized by a system function h(t). 
then the output is given by x(t) and is commonly symbolized as in Fig. 1. 





Figure 1: Input-output representation of a system in the time domain. 


For many reasons, one often transforms to the frequency domain (Fourier domain). Historically 
the most important reason has been that in that domain the relationship between input, output, and 
system function is simple for time-invariant systems In addition, and perhaps more importantly, 
the reason for going into the Fourier domain is that one can gain insight into the nature of the 
solution. If the input, system and output transforms of the time functions are given by F(uj) t H (a;) 
and A'(o;) respectively, then the input-output relations are symbolized as in Fig. 2. 

However, we have argued and have shown that this standard formulation does not fully describe 
what is happening if the system or signals are nonstationary, and that a significant simplification 
occurs if we have input-output relations in the time-frequency plane, as symbolized in Fig. 3. 
where Cf(t,uj) and C x (t,uj) are the input and output time-frequency distributions. The box with 
the question mark is meant to symbolize the time-frequency system function; we believe we have 
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F(oj) 


H(co) 


-► 



Figure 2: Input-output representation of a system in the frequency domain. 
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Figure 3: Input-output representation of a system in the time-frequency domain. 


solved the system approach as exemplified by the question mark, and have solved a number of 
problems to show the effectiveness of the method. In addition we point out that this approach 
lends itself to powerful approximation methods that give considerable insight into the nature of 
time-variant systems. 

Transformation of input-output relations into phase space. We write the input-output relations 
in the form 

L[t\x(t) = f{t) ( 1 ) 

where L is a linear (possible time-dependent channel). The system function can be a random 
function as we will describe consequently. For the governing system we take the linear time-variant 
system 


. . d n x(t) . , 

an (0 dt n + a "-l(0 dtn -1 ■■■ + ai(t) 

dt +a 0 {t)x{t) = f(t) 

(2) 

, d H , x tf 1 " 1 

d 


L(D,t) — a n(t) dtn +a n _i(*)^ n _, 

■ + a\ (£) — -1- ao(0 

(3) 

b 

II 

Sri ft. 


(4) 


For the moment, for the time-frequency distribution we take the Wigner distribution defined by 

Wi^c(t,u) = T Jx*{t - \r)x(t + ±T)e~ lT “dT (5) 
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We now state our main result: We write Eq. (2) in the following form 


L(D, t)x(t) = f(t) 


then the input output time-frequency relations are 


(G) 


where 


L*(A, £)L{B,J r )W x>x (t,uj) = W f j(t^) 

(7) 

a 1 d 

A = -- 

2 dt 

Id 

D — —— + JLJ 

2 Ot J 

(8) 

„ 1 d 

£ ~2id^ +t 

„ 1 d 

( ( J) 


The above are the basis of our approach, however, we emphasize that we have done this in 
more general ways by developing input-output relations when partial differential equations are the 
governing equations, and also we have considered the random case. We now describe some of 
the main results we have obtained. For ease of readership we have tried to make each section 
independent of each other and have, hence, repeated some of the equations rather then referring 
the reader to other sections. 


2 Results 

2.1 Quasi-Stationary channels 

Going from the time-invariant channels case to the time-variant case is a big jump, and one naturally 
may consider the in between case, which we call “quasi stationary” or “locally stationary”. By 
“locally stationary” we mean a system that is not stationary, but which is close to stationary for 
an interval of time. We have derived a criterion for local stationary in terms of time-frequency 
concepts. The basis of our derivation is simple and physical, namely that if around a time point the 
properties of the distribution does not change for that particular signal, then we have a situation 
that is indeed stationary around that time. We have obtained the following criterion: we consider 
a time point to and if 

B 7 (£o,u;) ~ Wt(to)W UJ {(jj) in the neighborhood of to (10) 

then we say that the process is locally stationary around the time point to- 
Also, if in an interval around a particular frequency point ojo 

\V(t,u) 0 ) ~ Wt{t)W 0 ) in the neighborhood of cjq (11) 
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then we will say t hat the process is locally stationary in the frequency variable. One can also define 
local stationarity in a time-frequency region. If at a time frequency point £o,u;o we have that 

W (£q,^ 0 ) ~ Wt{t neighborhood of£oM) (12) 

then one would say that the process is stationary in a region. For this to be so, both Eq. (10) and 
(11) must hold. In addition we have defined a measure of local stationarity by how much the joint 
distribution deviates from a product form. We use 

e 2 = [ f° + (F(t 0 ,w) -W{t 0 )W{u)) 2 dtdu (13) 

J Jto—S 

and similar criteria can be written for other types of time-frequency distribution. We have studied 
model systems and have found the above criterion to be effective. 


2.2 Other Joint Representations 


While we have used the Wigner representation, it is important to also do it for other representations, 
in order to ascertain which distributions are appropriate for different situations. We have obtained 
results for a number of special cases that are of particular interest. 

Ambiguity function domain. The ambiguity function is defined by 

■4»(t,u/) =-^ J x*(t — r/2)x(t + r/2)e ldi dt (14) 

and has been the main tool for radar and sonar systems. We have shown that 


L* (A a , £ a )L(B a , jT a )A x (t,u>) = A f (t,u) 


where 


A "~ 2 i0 dr' 

S = 11 - 1 , 

” ide 2 ' 


^ 15 1 

— T—J + “T 

1 06 2 


(15) 


(16) 

(17) 


Two time auto correlation function equation. We have obtained the equation for the two-time 
autocorrelation function. The two-time autocorrelation is the simplest quantity that characterizes 
how a random process x(£) is correlated with itself at two different times t \, £ 2 * It is defined bv 

Rx{ti,t 2 ) = E[x{ti)x*(t 2 )} (18) 

where E[ ] signifies ensemble averaging. For the differential equations given by Eq. (2), where 
now x(t) is a random variable and f(t) is a random driving process, the autocorrelation function 
R x (t\, ti) satisfies 

L L* Rx(ti,t2) — Rf{t\,t2) (19) 
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Continuous Wavelet Transform (CWT). The CWT of a signal x(t) is defined as 



(20) 


where a > 0, —oo < b < Too. The equation for the CWT is 


L{A w )C x (a, b) = C f (a, b) 


( 21 ) 


where 



( 22 ) 


2.3 Variance of a time-variant process 


In studying time-variant channels it is important to understand by how much they are varying 


as a function of time. We have shown that the exact instantaneous variance of a nonstationary 
random process can be obtained from the Wigner spectrum representation. What is advantageous 
in the evaluation of the variance with the time-frequency approach, is that we are integrating 
a smooth function, that is, the Wigner spectrum. Also, time-frequency distributions are very 
efficient in concentrating the information of signals in well localized regions of the time-frequency 
domain. The combination of smoothness and the concentration property is highly advantageous. 
The variance, of a the process v(t), is 


°l (0 = E KW] 


(23) 


where we have assumed that E [v(£)] = 0 for all times. We now show that indeed one can use 1 the 
Wigner distribution to calculate it. In particular 



(24) 


and hence, 



(25) 


To illustrate we consider the specific example 


^ +0v(t) = f(t ) 


(26) 


We have previously obtained the following exact solution 



1 A r o _ 2 0 t sin2 cjt No 
7T 2/3 e LO 27T 


2n (3 2 + a r 


1 
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Figure 4: Estimated and theoretical variance of v(t) of Brownian motion. 


and using this one gets 


(T 


2 

v 



1 No _ 2 /ft s in2u;£ No 1 
7r 2/3 lj 2n /3 2 -f u/ 2 


No e~ m 
2n (3 2 4- cj 2 


(cos 2u;£ — o;//3 sin 2a;£) 


do; 


(28) 


which simplifies to 


_2 _ M) 

v(<) _ 2/3 


1 — e 


-20t 


t > 0 


(29) 


As an numerical example to confirm that the instantaneous variance a^(t) can be obtained 
from the Wigner spectrum we performed a numerical simulation. We have generated the random 
differential equation with 10000 realizations and computed the variance using the Wigner spectrum. 
In Fig. 6 we show the estimated variance and the theoretical one. The agreement is excellent 


2.4 Transients 

If a time-invariant system, possibly random, is turned on, it will have transients until it achieves 
steady-state. Often, the transients are very important. We have been able to prove two important 
properties regarding such systems: 

a) We have shown that if wc turn on the system at any finite time £cn then the spectrum at 
infinity t = oo equals the standard power spectrum that is obtained with standard techniques. 

b) Furthermore, we have been able to show that the time-frequency output can always be written 
as the sum of a transient part and a stationary part which corresponds to the power spectrum. 
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In general the equation for the Wigner spectrum can be written as 


d 2n 


d 2n- 


d 


J2ir 


~ + l>2n-l o. 2n _l + • • • + ^>1 ^7 + ^ 


W x (t,u) = u(t)^ 


'dt 2n -' n_1 <9f 2n - 1 . " l dt 

where u(t) is the step function. We take the Laplace transform of Eq. (30) to obtain 

[&2u<S~ n + &2n-l5 2n 1 

Now, it follows that 


... + b\S + &o] W x (s, uj) — — - 

J 1 7T S 


SM = lim W x (t,u>) = \imsW x (s,u) = ^ ^ \H{u )\ 2 

t-*OQ s —►() Z7T f>0 Z7T 


(30) 


(31) 


(32) 


where H (u;) is the transfer function. This shows that when t —* oo the Wigner spectrum approaches 
the stationary solution. Now, it is always possible to write the Wigner spectrum fT x (£,u;) as a sum 
of a stationary part, constant with time, plus a transient spectrum. We write 


Wx{t,u) = W s {t,ij) + W T (t,u) 

To show this fact, one decomposes IT x (.s,i»;) into a sum of partial fractions 
No. 7V 0 A & 


2n 


2*N(s) 


(33) 


(34) 


fl (b 2 nS 2n + b 2n-l« 2n 1 + ... + b\s + b 0 ) 2tt S b2 n S 2n + &2„-iS 2rt_1 + .. . + b\S + bo 

where A (= ^ = |i/(u;)| 2 ) is a constant, and N(s) is a polynomial in the complex variable s. One 
then may prove that the stationary and transient part of the Wigner spectrum are 


777 / \ No A Nq I 2 1 

W s(s,w) = —- = — \H(u>)\ - 
27r S ITT s 

Ws(t,u) = ^u(t)\H(u)\ 


Wx{s f u) — 


2 

%N(s) 


hnS 2n + b<2n-\S 2n 1 + ... + b\s + 6 0 


(35) 

(36) 

(37) 


In the time domain, 

W s (t,u) = ^u(t)\H(u)\ 2 (38) 

which is exactly the stationary spectrum. 

2.5 Extension to MIMO 

Very often systems have multiple inputs and multiple outputs (MIMO systems) that are governed 
by differential equations of the form 

jX(t) = M(t)X(t) + F (t) (39) 
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where now X(£) and F (t) are column vectors with n elements 



( X,(i) ^ 


( fi(<) ^ 


X(<) = 

X 2 (f) 

. m = 

m 

(40) 


y x n {t) j 


y fn(^) J 



and M{t) is consequently an n-by-n square matrix of complex coefficients. The derivative in Eq. 
(39) implies that every element of X(£) has to be differentiated with respect to time. Now we define 
the autocorrelation matrix of the solution vector X(£) as 


flx(*l,*2)=£[x(*l)X(*2) f 


(41) 


The autocorrelation of the input vector F (t) is defined in the identical way. We have shown that 
the equation for the autocorrelation of X(£) defined by Eq. (39) is 

a n o a 

- ^-RxM'(t 2 ) - —M(t 1 )R x + MMRxM'fo) = Rf (42) 

where we have written Rx instead of Rx(t 1 ,^ 2 )? to keep the notation short. 


2.6 Input-output relations for partial differential equations 


We have developed input-output equations for systems governed by partial differential equations. 
We do not give here the equation, they can be found in our papers. We argue here that our 
method gives insight into the fundamental nature of certain partial differential equations. Consider 
for example two famous differential equations, the Schrodinger free particle equation and the heat 
equation. 

Shrodinger free particle; a = h/(2m) (43) 

Diffusion equation; D — diffusion coefficient (44) 


dip . d 2 ip 
~&t = la dx* 
du d 2 u 

dt dx 2 


Much has been written as to the analogy of the two equations, and 011 the fact that one can be made 
equivalent to the other by considering an imaginary diffusion coefficient. We have transformed the' 
equations into phase space where both become real. In this case the phase space is position and k 
space. The results are, respectively 


dWj, 

dt 

dWu 

dt 


= —2 ka 


d 2 W, h 


dx 2 


D d' 2 \\\ 
~2 dx 2 


- 2 Dk 2 W u 


(45) 


(46) 
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where 

W u {x. k; t) = 2- J 

While the only difference between the two original equations. Eq. (43) and (44), is an /, the 
difference in phase space can be readily seen and studied, and indeed each of the terms has a 
physical interpretation. We have transformed other partial differential equations, and in each case 
the physical nature of the original equations is readily seen. Moreover, we have shown that solving 
partial differential equations in phase space is simpler. 


A 

2' 1 


A 

x + 2’ t 


- jkx d\ 


(47) 


2.7 A 11011 -white noise model 


We have developed an example for the non-white noise case that can be done analytically and can 
be used as a model example. Consider the random differential equation 

d 2 x(t) 0 dx(t) 


dt 2 +l3 dt t( i ) 


where £(£) is the noise operator that satisfies 


<me(0]+) = ^ f° 


(48) 


(49) 


It is called the quasi-classical Langevin equation and is standard in quantum optics. In such a case 
the autocorrelation function is 


R(t) = — [°° ue MT coth -^—duj 
K ’ 2tt _ 2 kT 


We rewrite it as 


with 


p(t) + /3p(t) - £(t) 




i? c (r) - 2 DZ / u/e iu,r coth ZuxLj 


h 

2 kT 


Using the standard Wiener-Khinchen theorem the power spectrum is given by 

S^(u) — 2D Zuj coth Zuj 

We note that as Z —► 0 

|ini R^t) = 2 DS{t) 


(50) 

(51) 

(52) 

(53) 

(54) 


which is the white noise case. Using our general approach we have that the differential equation 
for the Wigner spectrum is 


Id 2 a d 2 2 

4 W +P di +f3 + UJ 


_ [)% 

W p (t, u) ~ - uj coth Zuj 

7 r 


and we have found the exact solution 


W p (Llj) = 


1 DZ 


(3 2 + uj 2 7 r 


u: coth Zuj 


_2x t 

1 — e rn COS 2 ujt 


( 55 ) 


(56) 
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2.8 Dispersive channels with damping 


Channels are generally dispersive and also have damping. The fundamental nature of dispersive 
wave propagation is that different frequencies travel at different velocities, and hence one would 
expect that a transformation into a joint posit ion-wavenumber representation would be well suited 
to study such dispersive propagation. Using our method we have shown that a pulse can be 
propagated in a very simple way without solving the governing differential equation. The classical 
method for solving a wave equation with constant coefficients, 

N ™ M r) m n 

(57) 


h 

«=0 m=0 


dx 11 


is to substitute e lkx lLJt into it, to obtain a relation between k and u 

N M 

yi On(-w) w = b m {ik) m (58) 

n=0 rn=0 

One then solves for uj in terms of k to obtain the dispersion relation uj = uj(k). There may be more 
than one solution and each solution is called a mode. The general solution for a particular mode is 


u(x 


with 


,o= ;/g/ s(M)e ' 


J u(x, 


kx—iu/(k)t dj ^ 


S(k, 0) 


0) e~ ikx dx 


\p2/n 


If one defines the time dependent spectrum by 

S{k,t) = S{k, 0)e~ iuj{k)t 

then u(x,t) and S(k,t) form Fourier transform pairs between x and k 


u(x,t) = 


wJ s{k ' t)e " 


dk 


V2n . 

S(k, t) = -L= [ u{x , t) e~ lkx dx 

v 27r J 

For a complex dispersion relation 


uj{h) — ujfi(k) + iujj{k ) 


the group velocity, v g {k ), is defined by 


v g {k) = J R {k) 


(59) 

(60) 

(61) 

(62) 

(63) 


(64) 


(65) 
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The posit ion-wave number Wigner distribution at time t is defined as 

W(x,k, t) = J u*(x — ^r,t)u(x + \r,t) e~ lTk dr (6G) 

where u(x , 0) is the input. We obtained an approximation to \V(x, fc, t) in terms of W(x, k , 0) 

W a fa k, t) = e 2u,/(fc) * W(x - (fcjt, fc, 0) 


where W a is the approximate Wigner distribution. This is a remarkably simple result and besides 
being easily applied, this approximation is interesting and insightful because it is physically inter¬ 
pretable. In particular, each local point in phase space moves at constant velocity given bv the 
group velocity, with attenuation governed by ujj(k). We have tested this result on many exactly 
solvable examples, and it works much better than current approximation methods, such as the 
stationary phase approximation. 

Ambiguity function. Extending this approximation approach to other representations is po¬ 
tentially useful. The reason is that while different representations have properties in common, 
sometimes different representations are more tractable. Of particular interest is the ambiguity 
function domain 

A/(0, r, t) = J u*(x — ^r,t)u(x 4* \r,t) e l0x dx (67) 

We note that the ambiguity function and the Wigner distribution are a Fourier transform pair. We 
have been able to show that 


M(e,T,t) = f^J J A/(0 1 r',O)e ! ' fc(r - T ' ) e ^*(H-|0)-u,(fc-±0)] ^ dk (68) 

This is exact. We have also been able to show that the approximate ambiguity function is 

M a (0,T,t) = f- J J e^ ,(k)t M{6, t',Q) e ik(r - T,) e' tv * {k)0 dr'dk (69) 

Spectrogram. The spectrogram is perhaps the most widely used distribution and goes under 
many names, for example in acoustics it is called the u lofargram”. It is often used to refer to 
the time versus frequency plot of signal energy obtained by a filter-bank analysis. Since we are 
considering the phase space of position and wave number, we define the “local” (or “short-space”) 
wave number transform by 

~ ikx ' dx' 


S x {k,t)= \ _ [ u{x\ t)h*(x — x')t 

v27t J 


(70) 


where h(x) is a spatial window function that is typically narrower in extent than the wave u(.i\ t). 
peaks at x = 0 and tapers to zero for \x\ > 0, although this is not required in what follows. The 
position-wave number spectrogram is then given by 

1 


P sp (x,k,t ,) = | S x {k,t) | = 


v/2t r 


u(x\ t)h*(x — x')e lkx dx 


(71) 
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We have found the following approximate solution 

P sp (x, A, t) * T J e 2u ' , (( fc,+k ")/ 2 ) < S(A',0)5*(A", 0) /T(A - A:') 

//(A - A") c -«*-«s»((*'+*")/ 2 )0(* w -*') dfc'dfc" (72) 

Notice now that this no longer can be put in the spectrogram form, that is, as a magnitude- 
square of some integral function. 


2.9 Position-Wavenumber Approximation for Random Dispersive Channels 


We derived an approximation of the position-wavenumber (x-A:) Wigner spectrum of the random 
channel, and its impact on a propagating wave. When there is no dispersion and the channel is 
deterministic, the pulse does not spread as it propagates However, in a random channel, the pulse 
spreads in the mean even in the absence of dispersion. Let u(x, 0) be an initial deterministic pulse, 
propagating in a random channel with Wigner spectrum W / t (x, A;; t). The propagating wave u(x, t) 
is therefore also random, and its Wigner spectrum W u (x,k; t) is given by [1, 7] 


W u (x, k; t) 


J W u {x - x', A; 0) W h (x\ A; t) dx' 


(73) 

(74) 


where W u {x, k: 0) is the Wigner distribution of it(x,0) and P{.} denotes the expected value. The 
approach is to extend the deterministic case to the random case by introducing random parameters, 
and then ensemble averaging the deterministic approximation. 

As a case/example, we consider a dispersive channel with random speed c and exponential 
attenuation parameterized by the variable b , with joint probability distribution P(6, c). For this 
case, the approximate Wigner spectrum of the channel is given by, 


J e- m p db 

where v g (k) = cv\(k) is the group velocity, parameterized by the random sound speed c and the 
deterministic function v\ (fc). The dispersionless case corresponds to v\ (k) = 1. 

Let the joint probability distribution function P{b, c) be independent, with sound speed r de¬ 
scribed by arbitrary distribution P c {c) and the damping parameter b given by an Erlang distribut ion. 


Pb(b) = 


1 


A n b n - l e~ x \ 6 > 0 


(n- 1)! 

The Wigner spectrum of the channel is then approximately 

W h (x,k-,t) ( A v * 1 


\ + ktJ |£ uj (A) | \tv\(k) 


(75) 


(76) 
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To obtain the approximate Wigner spectrum of the propagating wave, we convolved this channel 
spectrum with the Wigner spectrum of the initial wave. Note that as the wave evolves the channel 
increasingly attenuates the propagating wave, and the statistics of the wave in x are distributed 
according to Pc (r^riy)- If we take P c (e) to be Gaussian with mean cq and variance a 1 , then the 
approximation is 


W fl (x.k; t) 


X + htJ ^/2irt 2 vf(k) a 2 


(x 

2 t* vf(k) a 2 


(77) 


Note that at / — 0 the wave has not propagated, and the channel approximation becomes <S(a;), 
which is a satisfying result in that this result is exactly satisfied by the approximation at t = 0. 
This is not true of other standard approximations, such as the stationary phase approximation, 
which is not accurate for small times. Also, note that even when there is no dispersion {vi(k) = 1). 
the channel still induces spreading on average in the wave, unlike the deterministic case for which 
the channel Wigner spectrum becomes 5(x — c£) in the absence of dispersion (and damping). 


2.10 Invariant Features for Classification 

As a pulse or wave u(x , t) propagates, it can change its shape due to factors such as frequency- 
dependent attenuation and dispersion by the propagation channel. This can have a deleterious 
impact on detection and classification. Hence one asks is there are features that are invariant to 
propagation effects. We have derived new central-like moments that are invariant to dispersive 
propagation as we know explain. 

Let u(0, t) denote the initial pulse as a function of time t , generated at position x — 0. Then, 
the pulse at a subsequent position is given by [4, 5, 9, 10] 

u(x,t) = — L= f F(0,u;) e~* u,t duj = L_ f F{x,uj) e~^ 1 doj (78) 

v 27 T J v 27 t J 

per mode, where F(0,u;) is the spectrum of the initial pulse, 

F(0,w) = —L= [ u{0.t) e>“ l dt (79) 

\/27r J 

K(oj) is the dispersion relation expressed in terms of A: as a function of uj. If the dispersion 
relation is complex, then there is damping, resulting in frequency-dependent attenuation as the 
wave propagates. Also, we define the amplitude and phase of the wave and its Fourier transform 
by 

u(x, t) = A(pc, t)e^^ x ^ 

F(x, w) = D{x,uj)e j ^ x ^ 
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(80) 

(81) 












Note that 


F(0,u) =B(0.w)^ (m 


(82) 


and therefore 


F{x,u) = F{ 0,u)e jK{uj)x 

= D(Q,cj)e j ' i ’ {0 ' u,) e jh ' Mx 


(83) 

(84) 


It follows that, for the real dispersion relation K(u>) 


B(x,lj) — B( 0,cj) 

^{x.uo) = ^(0, lj) -f K(u)x 


(85) 

( 86 ) 


The eentral-like moments that are invariant to dispersive propagation are given by 



(87) 


where tg{x,u) = — is the group delay of the wave. These moments are similar to central 

temporal moments, but with an important difference. To see the similarity, consider ordinary 
central temporal moments of order 7i, given by 



( 88 ) 


(89) 


where Eq. (89) follows from Eq, (88) by substituting in for u(x, t) in terms of its Fourier transform. 


and (t) is the mean time. (Note that for n = 2 this is the duration moment.) We see that the 


moments A n (x) are similar to ordinary central temporal moments, except that rather than being 
centered about the average time ( t ) x , they are centered about the group delay t g (x,u) of the wave. 

This centering about the group delay affords these moments with a property that is particularly 
attractive to their use as potential features for classification, namely, they do not ehange with 
propagation, even though the wave it(x, t) does. That is, they are independent of location x, [2, 6], 


A n {x) = A n (0) 


(90) 


for real dispersion relations. Thus, if differences are observed in these moments for two differ¬ 


ent received pulses, the differences reflect differences in the source rather than differences due to 


propagation (i.e. dispersion) effects. The moments A n (x) ean be equivalently calculated in the 
time domain [2, 6], which is beneficial for computational purposes. In particular, the equivalent 
time-domain formulation is 



(91) 
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Also, as defined, the odd-order moments are identically zero. Hence we use a one-sided integral 
over t in Eq. (91) to obtain non-zero odd-order dispersion-invariant moments. 


It is important to appreciate that the results above are for a mode. Since in general each 
mode travels with a different group velocity (i.e., each mode has a different dispersion relation). 


the spectral and A n moments of the total wave, which is comprised of the sum of modes, change 
with distance, even though the spectral and A n moments of each mode are invariant to dispersive 
propagation. Thus, inode separation prior to calculating the moments is important in order to 
obtain invariant moments. However, even without mode separation, the moments A n will be less- 
variable with propagation distance than will ordinary temporal moments, and may therefore still 
be better features for classification. 

We have compared classification performance to that of ordinary temporal moments on nu¬ 
merical models of acoustic scattering from steel shells in a Pekeris waveguide. We evaluated the 
two-class problem of distinguishing a sphere from a cylinder, and the more challenging problem 
of distinguishing between two different cylinders, via their backscattered echoes at various prop¬ 
agation distances. The classification utility of our moment features was assessed by computing 
receiver operating characteristic (ROC) curves for each two-class problem. We also considered an 
approach based on a correlator receiver. We applied our dispersion-invariant approach to define a 
dispersion-invariant correlation coefficient for use as a classification feature, and have shown the 
detrimental effects of dispersion on the ordinary correlation coefficient. We calculated ROC curves 
and have shown the significant performance gains over the ordinary correlation coefficient for the 
two-class simulations. 

2.11 Fields 

We developed a direct method using phase space distributions, such as the Wigner distribution, to 
study wave propagation in a dispersive medium. We derived an explicit evolution equation for the 
Wigner distribution and methods to solve it. This has been described in the previous sections. We 
have also developed the method to the evolution of noise fields in a dispersive medium. For linear 
wave equations with constant coefficients the solution is given by 



(92) 


where S(k. 0) is the initial spatial spectrum which is obtained from the initial pulse by way of 



(93) 


If we define 


S{k,t) = S{k, Q)e~ lu{k)t 


(94) 
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iii which case 

u{x,t) = -L [ S{k,t)e ikx dk 

V 27T J 

S(k,t) = -— / u(x, i) e _l ^ x rfx 

v 27r J 

Differentiating Eq. (92) with respect to time and multiplying by i, gives 

i^u{x:t) = [ u{k)S{k, 0)e~ iu{k)t e ikx dk 

ot v27T J 


(95) 

(90) 

(97) 


Now, in general, the Fourier transform of the product of two functions can simplified bv the use of 
the following relation, 


/ wsW" <*-/' (\£) ^ dk - i Gl) / **> 


e ikx dk 


Applying this to Eq. (97) we have that 

.d 


We rewrite this as 


i—uix, t) = uj 
dt. ’ 


.8 


id_ 

i dx 


u(x , t) 


i—u(x,t) = u(IC)u(x,t) 


where /C is the “wave number operator” defined by 

We derived the equation of motion for the Wigner distribution 


ij t W(x,k,t) = 


uj[k+ hiL) ~ u ( k ~h^ 


W{X, k, t) 


(98) 


(99) 


( 100 ) 


( 101 ) 


( 102 ) 


We note that this equation is derived under the assumption that uj is real. In Eq. (102) what is 
meant by c o{k + is that we substitute the operator k + in u(k ). 

Noise fields : The problem we have done is the following. We consider the case where at a 
spatial point, xq, a noise field is produced as a function of time and we ask what the field is at a 
later time. One has that. 


u(x. 






(103) 


where K(u>) is the dispersion relation as a function of uj and where F(uj) is the initial time spectrum 
at x = a*o, 

e iu,l dt (104) 
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By defining 


F(u,x) = F(u,x 0 )e lK ^ )x 

where 

F(u,xo) = F(u>) 

we have that 

u(x,t) — j_ f F(uj, x) e~ lujt (Lo 
v27r J 

F(lj,x) = \— [ n(x, t) e lu;t dt 

v 27r J 

In addition it is convenient to combine Eqs. (103) and (104) 


(105) 

(106) 


(107) 

(108) 


u(z, t) 





(109) 


Let us now study the statistical properties of u(x, t) considered as a random field. We first ask the 
following question. Suppose the noise field is stationary at position of generate xo- will it be time 
invariant at other positions? In particular consider 


u*(x , t) u(x, t 


t) = jj u*(x 0 ,t')u(x 0 ,t") 

c iK(«y)(x-i 0 dt'cLodfdu 


e -1K (a?) (x - x 0 ) +iu>t e - iujt' 


(no) 


and take the ensemble average 

E[u*(x,t)u(x,t + r} = ^2-^ J J Elu^xo,t')u(x 0 X)}e~ iKM{x ~ xo)+iut e- iu,t ' 

e iK(«,')(x-x 0 dt'dwelt" (ko' ( 111 ) 

where we have assumed that K(u>) is a deterministic function. Suppose we define R by 

R(xQ,t” - t') = E[u*(xq, t , )u(xo,t ,, )\ (112) 

where we have explicitly assumed stationarity. Hence 

£[u*(i,(Mi,f + r)] = ( 2) j R(xq, t" - i , )e‘ iJC (")(*"* 0 ) +<w ‘e-"*' 

e iK(u,')(x-x 0 )-xu\t+r) e iu>'t'' dt > d ^ dt n d j (H3) 

We have been able to show that this reduces to R(x o,£ /; ). Thus we see that if the process is 
stationary at position xq it will be stationary at any point in space. We now consider the same 
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issue from a different perspective. We define the Wigner distribution in t,u> of u(x, t) at a particular 
position x by: 

W x) = — [u* (x,t — r/2) u (x, t + r/2) e~ lurr dr (114) 

2tt J 

One defines the Wigner spectrum by 

W(t,u\x) = [ E[u(x, t — t/2)u(x, t -f- r/2)]e“ lTu; dr (115) 

2n J 

We have shown that 

W(t,u; x) = 2- J J W (t\^x 0 )e-' e{t '- t) e- l{x - xo)[k ‘^ +d/2) - k ^- d/2)] dtUW (116) 

which relates the Wigner distribution at position x to the Wigner distribution of the initial pulse 
at x — £o. We now take the ensemble average of both sides of Eq (11G). 

W x) = 2- J J x 0 ) e -*«(«'-O e -i(x-x°)[fc*(w+e/ 2 )-fc(w-e/ 2 )]^ d6 , (117 ) 

Now, let us assume that W xo) is independent of time which is the condition of a stationary 
process. Let us call it W (a;; x*o) . We have 

W(t,w; x) = j I W(oj; 0)S(e)e- l ^ t '- t) e- lx ^ {uJ+9/2) - k ^- d/2)] dt'de = W (u: 0) (118) 

which shows that the Wigner distribution is also stationary at position x which is consistent with 
what we derived above. 

In addition we consider the same problem using the approximate Wigner distribution. We have 
shown that an very accurate approximation can be obtained by way of 

W (£,u>; x) & W (t - t 9 (cj) (x - x 0 ),u;;xo) (119) 

where r g (<jj) is the derivative of the dispersion 

r g (u) = k'( u») = ^ (120) 

and called the “group slowness”, or “unit transit time”. 

Now, take the ensemble average of Eq. (118) 

W(t,uj;x) ^ W(t - T g (a;) (x - x 0 ),cj;xo) (121) 

and it is clear that the approximation also satisfies Eq. (118). We emphasize that the above 
considerations is for the case where the dispersion relation is deterministic. In many cases of course 
that is not the case, and such situation is under development. 
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